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(57) ABSTRACT 

This invention is drawn to an autonomous navigation system 
using Global Positioning System (GPS) and magnetometers 
for low Earth orbit satellites. As a magnetometer is reliable 
and always provides information on spacecraft attitude, rate, 
and orbit, the magnetometer-GPS configuration solves GPS 
initialization problem, decreasing the convergence time for 
navigation estimate and improving the overall accuracy. 
Eventually the magnetometer-GPS configuration enables the 
system to avoid costly and inherently less reliable gyro for 
rate estimation. Being autonomous, this invention would 
provide for black-box spacecraft navigation, producing 
attitude, orbit, and rate estimates without any ground input 
with high accuracy and reliability. 

24 Claims, 7 Drawing Sheets 
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AUTONOMOUS NAVIGATION SYSTEM 
BASED ON GPS AND MAGNETOMETER 
DATA 

CROSS-REFERENCES TO RELATED 5 

APPLICATIONS 

This application now formalizes and incorporates herein 
by reference Provisional Application Serial No. 60/301,363, 
“Autonomous Navigation System Based on GPS and Mag- 
netometer Data” Richard R. Harman, et al., filed on Jun. 25, 10 
2001. Applicant claims the priority date thereof under 35 
U.S.C. 119(e). 

ORIGIN OF THE INVENTION 

The invention described herein was made by employees 
of the United States Government. The invention may be 
manufactured and used by or for the governmental purposes 
without the payment of royalties thereon or therefor. 

FIELD OF THE INVENTION 20 

The present invention is directed to a method and appa- 
ratus for an autonomous navigation system for estimating a 
set of satellite navigation parameters based on measure- 
ments from Global Positioning System (GPS) and magne- 25 
tometers. 

BACKGROUND OF THE INVENTION 

GPS has become a standard method for low cost onboard 
satellite orbit determination. The use of a GPS receiver as an 30 
attitude and rate sensor has also been developed in the recent 
past. GPS is emerging as the favored technique that provides 
sufficient accuracy for most requirements. 

However, GPS requires expensive gyroscopes (gyros) or 
star trackers to be fully functional and provide accurate 35 
attitude data. In addition, GPS is subject to initialization 
delays and signal drop-outs that require back up systems to 
fill in during the loss of GPS data. 

SUMMARY OF THE INVENTION 40 

This invention presents a method and apparatus for pro- 
cessing GPS phase, pseudo -range, and range rate data as 
well as magnetometer data to produce estimates of attitude, 
angular rate, position, and velocity for a spacecraft, espe- 
cially a low earth orbit (LEO) satellite. 

Magnetometers are a low cost and reliable technology that 
have been used to determine satellite orbit and attitude based 
on the earth’s magnetic field. Magnetometers always output 
data as they are not subject to occultation or tracking 5Q 
problems. 

Although GPS is relatively a low-cost sensor, providing 
accurate orbit estimates, with coarse attitude estimates, GPS 
based systems are subject to initialization problems, and 
must rely on costly systems such as gyros. However, GPS 55 
data is available continuously throughout the orbit and can 
produce more accurate orbit, attitude, and rate estimates by 
combining it with magnetometer data. 

The present invention solves the drawbacks associated 
with GPS -based systems as mentioned above, by coupling 60 
GPS with reliable and low cost magnetometers. By process- 
ing the orbit and attitude solutions in the same algorithm in 
a single self-contained unit, weight, size, and power con- 
sumption are reduced. Completely autonomous satellite 
navigation is thus possible with reduced weight and cost. 65 

The method used is an extended Kalman filter (EKF), 
blended with a pseudo-linear Kalman filter algorithm, devel- 


2 

oped to provide estimates of attitude, orbit, and rate esti- 
mates using magnetometer and GPS data. The magnetom- 
eter based EKF can converge from large initial errors in 
orbit, attitude, and rate. Combining the magnetometer and 
GPS data into a single EKF provides a more robust and 
accurate system that produces estimates of attitude, angular 
rate, position, and velocity for a spacecraft. 

Some publications listed, relating to the present invention, 
are incorporated by reference: 

1. Deutschmann, J., R. Harman, and I. Bar-Itzhack, “An 
Innovative Method for Low Cost, Autonomous Naviga- 
tion for Low Earth Orbit Satellites”, Proceedings of the 
AAS/AIAA Astrodynamics Specialists Conference , 
Boston, Mass., Aug. 10-12, 1998. 

2. Oshman, Y, and F. L. Markley, “Spacecraft Attitude/Rate 
Estimation Using Vector — Aided GPS Observations”, 
IEEE Transactions on Aerospace and Electronic Systems, 
Vol. 35, No. 3, July 1999. 

3. Deutschmann, J., and I. Bar-Itzhack, “Comprehensive 
Evaluation of Attitude and Orbit Estimation Using Real 
Earth Magnetic Field Data”, Proceedings of the 11 th 
Annual AIAA/USU Conference on Small Satellites, 
Logan, Utah, Sep. 15-18, 1997. 

4. Bar-Itzhack, I., “Classification of Algorithms for Angular 
Velocity Estimation”, Journal of Guidance, Dynamics, 
and Control, Vol. 24, No. 2, March-April 2001. 

5. Deutschmann, J., R. Harman, and I. Bar-Itzhack, “A LEO 
Satellite Navigation Algorithm Based on GPS and Mag- 
netometer Data”, Proceedings of the CNES 15 th Interna- 
tional Symposium on Space Flight Dynamics, Biarritz, 
France, Jun. 26-30, 2000. 

6. Brown, Robert Grover, and Patrick Y. C. Hwang, Intro- 
duction to Random Signals and Applied Kalman 
Filtering, (3 rd ed.), John Wiley & Sons, 1997. 

7. Wertz, James. R, editor, Spacecraft Attitude Determina- 
tion and Control, D. Reidel Publishing Company, 1984. 

8. Shorshi, G. and I. Bar-Itzhack, “Satellite Autonomous 
Navigation Based on Magnetic Field Measurements”, 
Journal of Guidance, Control, and Dynamics, Vol. 18, 
No. 4, July-August, 1995. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 shows a block diagram of an autonomous naviga- 
tion system based on GPS and magnetometer data. 

FIG. 2 shows the position errors about each axis of the 
Simulation Case 1. 

FIG. 3 shows the velocity errors about each axis of the 
Simulation Case 1. 

FIG. 4 shows the angular attitude error about each of the 
body axes of the Simulation Case 1. 

FIG. 5 shows the rate errors about each axis of the 
Simulation Case 1. 

FIG. 6 shows the result of applying the method of the 
present invention to RSS Angular Velocity Error During 
First 150 Minutes, Measurements from One GPS Satellite of 
the Simulation Case 2. 

FIG. 7 shows the results of applying the method of the 
present invention RSS Angular Velocity Error During First 
150 Minutes, Measurements from Two GPS Satellites of the 
Simulation Case 2. 

FIG. 8 shows the results of applying the method of the 
present invention Estimated Angular Velocity About Y Axis 
(Estimation Approach) of the Simulation Case 2. 

FIG. 9 shows the results of applying the method of the 
present invention Estimated Angular Velocity About Y Axis 
(Derivative Approach) of the Simulation Case 2. 
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FIG. 10 shows the results of applying the method of the 
present invention RSS Attitude Error During First 150 
Minutes, Measurements from One GPS of the Simulation 
Case 2. 

FIG. 11 shows the results of applying the method of the 
present invention RSS Attitude Error During First 150 
Minutes, Measurements from Two GPS of the Simulation 
Case 2. 

FIG. 12 shows the results of applying the method of the 
present invention RSS Position Error, Measurements from 
One GPS of the Simulation Case 2. 

FIG. 13 shows the results of applying the method of the 
present invention RSS Position Error, Measurements from 
Two GPS of the Simulation Case 2. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS OF THE INVENTION 

FIG. 1 shows a block diagram of an autonomous naviga- 
tion system 10 based on data from GPS and magnetometer 
26. The system includes GPS receiver 18, magnetometer 26, 
process electronics and interface, and the navigation algo- 
rithm 20 embedded on a process board containing naviga- 
tion algorithm 20 , and power supply 12 (either external or 
internal). 

The GPS/magnetometer navigation system is based on an 
extended Kalman filter algorithm, blended with a pseudo- 
linear Kalman filter. The general EKF equations are: 

(!) 

Z A+1 = A a+ i(JC(^+i)) + a A+1 (2) 

The general equations for the pseudo-linear Kalman filter 
are: 


4 

include j2 and other higher order terms, using equation (1). 
The clock errors are propagated according to a two state 
random-process model. The quaternion is propagated 
numerically using the kinematic equation. Finally, the rota- 
5 tion rate is propagated using Euler’s equation and b is 
propagated with the equation expressing the time rate of 
change of a vector with respect to a rotating coordinate 
system. The two equations are augmented as 


10 

(h 


/ 1 [(/£> + £)*] 0" 

' 01 

+ 



b 


B 0. 

_b_ 


U_ 


where 

15 B=[b m x]=cross product matrix containing the measured 
magnetic field vector, b m 
I=inertia matrix 

h=angular momentum of reaction wheels 
20 T=external torques on the spacecraft 

U=Db r 

25 D=direction cosine matrix which transforms from inertial 
to body coordinates 

b r =reference magnetic field vector in inertial coordinates, 
computed from the estimated position 
Equation (5) is treated as a ‘pseudo -linear’, equation in 
30 which the dynamics matrix contains the current best esti- 
mate of the rate, 

&>. 

35 

The filter covariance matrix is propagated according to 


m = F(X(t k ))*X(t) + m 


(3) 


P k+ 1 (-) = <*>k {x k (+)) A (+)<*>[ (£* (+)) + Qk 


( 6 ) 


£t+i - ^ife+i *£t+i +a Jt+i 


(4) 40 


The output filter state is given as 

x t =[r t , v T , c T , a T , w r , b T ] 

45 

where R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver clock 
errors, q is the attitude quaternion, cq is the rotation rate, b 
is the magnetic field in body coordinates. 

The filter processes data from the magnetometer, a 50 
2-baseline GPS receiver, and reaction wheels (used in Eul- 
er’s equations). The GPS data, as mentioned, consists of 
range, range -rate (if available, or reliable), and phase mea- 
surements. Nominally, the filter first propagates the state 
estimate and filter covariance and then performs an update 55 
with the magnetometer data followed by an update with the 
GPS data. Data from one, two, or more GPS satellite 
observations are included in each update. The propagation of 
the state estimate and the filter covariance will be described 
first. The description of the state and covariance update will 60 
follow, including the development of both the magnetometer 
and GPS measurement matrices. 

Filter Propagation 

The dynamics for the position and velocity, the 
quaternion, and the rotation rate and b are treated as 65 
uncoupled. The position and velocity are propagated 
numerically with a simple two-body model which can 


<D=the state transition matrix=e jFAr . 

F =fflfX 

Q=matrix of standard deviations of the white noise pro- 
cesses driving each of states. 

In general, the linearized error dynamics are given as 


x{t) = F*x(t ) 


(7) 


Where the dynamics matrix, F, is computed for the 
position, velocity, and quaternion. The error state is repre- 
sented by x for the position, velocity and quaternion. The 
rate and b state equation is linear and therefore the dynamics 
matrix is given in (5) (and is used to propagate the error in 
the actual states cq and b). The matrix F is partitioned as 


F orbit 0 0 0 

0 F clock 0 0 

0 OF attitude 0 

0 0 0 F rate , b 


where: 


forbid 


(8) 
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F clock — 


1 O' 
0 0 


F attitude ~ 

[ r l [{io± + h)x ] o 

B 0 


The innovations for the rate and b states is then 


~-b m -H ra 


i(-) 


( 15 ) 


5 

where 


K rat ehjc+i (-) = the estimated rate 


The filter uses the ‘multiplicative’ approach in the attitude io 
estimation, therefore the error state for the attitude is a 
vector of attitude errors. 

Filter Update 

The EKF state and covariance update equations are given 
as 

=i* +1 (-)+* 1+ ik +I -M-W-)] 

Pm (+) = [/- K k+l H m (t k+i (-)] 

The Kalman gain is computed according to 


and b components of the state vector, and the other param- 
eters are as defined above. This defines the estimation 
approach for the magnetometer update. 

GPS Update 

15 The measurement matrix for the position, velocity, and 
clock states (combined) is a standard model based on range 
and range -rate data (actual GPS measurements would be 
pseudo -ranges and delta-ranges, but for the purposes of 
testing, range and range-rate measurements were generated). 

20 The measurement matrix is obtained by linearizing the range 
and range -rate equations that are functions of the spacecraft 
position and velocity, respectively. (If the range -rate is not 
used, the matrix only contains the range development.) The 
resulting measurement matrices are 


K k+ 1 = 


( 11 ) 


25 




c x 

0 


Cy C z 0 

0 0 c* 


0 

c y 


0 1 0 
c z 0 1 


(16) 


where H=fh/fX 30 

The measurement matrix, H, is partitioned as H=[H orfetf , 

H attitude? H rat e,b\- The individual measurement matrices are 
defined according to the two measurements used, namely the 
magnetometer measurement and the GPS phase, range, and 
range-rate (if used) measurements. First, the development of 35 
the magnetometer measurement matrix will be presented, 
followed by the derivation of the measurement matrix for 
GPS. 

Magnetometer Update 

The measurement matrix for the orbit terms, H orfclY was 40 
developed in Reference 8 in terms of spherical coordinates. 
The additional relationship to express the partial derivative 
in terms of the position and velocity is straightforward. 
Obviously, the portion of Yi orbit , related to the clock errors is 
zero for the magnetometer. The measurement matrix for the 45 
attitude states is given as 

H attitud A-knA ( 12 ) 


Where c x , c y , c z , are the appropriate direction cosines 
between the spacecraft and the GPS satellite in view. The 
innovations are computed by differencing the measured 
range and range -rate with the computed values (computed 
range and range -rate are based on the filter estimated posi- 
tion and velocity), given as 

P pc (17) 

Z k+ 1 = 

Pk + 1 Pck + 1 

where p,pC=measured and computed range, respectively 
and p,pC=measured and computed range-rate, respectively. 

The measurement matrix for the attitude states is found by 
linearizing the relationship between the phase measurement 
and the line of sight vector to the GPS satellite in view of the 
antenna. Assuming a 2-baseline GPS receiver, the phase 
measurement for a given baseline can be written as 

y j r*j T Ds H w i (18) 


The innovations for both the orbit states and the attitude 50 
states is given as 

=tJ>b r (13) 

where 55 

D=(I-[ax])D=the direction cosine matrix computed with 
the filter estimated quaternion; and 
D=true direction cosine matrix which transforms from 
inertial to body coordinates. 60 

Finally, the measurement matrix, H rate b , for the rate and 
b is defined by the following linear relationship. 


[0 /] 


: Fl r at e jo 


( 14 ) 


65 


i=GPS satellite designation 
j=antenna designation (j=l, 2) 
y=the phase measurement 

a 7 =unit vector representing the baseline direction for 
antenna j 

s /r =unit vector to GPS satellite i, in inertial coordinates 

v~white sequence measurement noise 

The computed phase measurement is given as: 

y .. = ajfis . d9) 

-Jl J n 


The innovations, Z k+1 , is computed by subtracting (19) 
from (18). The equation can then be recast to produce the 
necessary linear relationship between the and the 

attitude error vector, defining the attitude measurement 
matrix for the phase measurements 
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where 


2 / ryji-yjr-afV^m^+vrHatSk+Vi 


( 20 ) 


*ai,m ~ “ 

U^-yl-yl \ 

measured unit vector to GPS satellite i in body coordinates; 

and a= attitude error state 

Finally, the measurement matrix for the rate states is 
developed by first differentiating the cosine of the angle 
between the baseline and the line of sight to a given GPS 
satellite gives (dropping the satellite and antenna 
designations) 

-sin(0)0 = a T D& r + a T bs r + a T D& r (21) 

where the quantities are as defined above. Substituting in the 
kinematic equation for the direction cosine matrix, and 
noting that 

a = 0 


gives 


-sin(0)$ = -a T [oix]Ds. r + q? D s r 


(22) 


30 


Defining the ‘measurement’ to be 
ya, = a T Ds + sin(0)$ 


(23) 


35 


8 

Augmentation of GPS Phase Measurement 

Now the estimation approach is applied to the GPS phase 
measurements, which contains either the phase measure- 
ment states for one GPS observation, two, or more. 

The output filter state is given as 


: [fl r , V T , C T , q T , (J, b T , z T ] 


(27) 


Filter Propagation 

Again the dynamics for the position and velocity, the 
quaternion, and co, b, and z are treated as uncoupled. The 
position and velocity are propagated numerically with a two 
body gravity model which includes the j2 term. The clock 
errors are propagated according to a two state random- 
process model. The quaternion is propagated numerically 
using the kinematics equation. 

The rotation rate is propagated using Euler’s equation. 
The states b and z are propagated with the equation express- 
ing the time rate of change of a vector with respect to a 
rotating coordinate system. The equation for the rate and b 
are presented above. 

Following is the development of the equation expressing 
the time rate of change of the elements of z. 

The expression for a single GPS phase measurement is 
given as 


Z ij~ a i' S j ( 28 ) 

where 

phase measurement; 

a —unit vector in direction of baseline i(i=l,2); and 

s 7 —unit vector in direction of a GPS satellite (j=l or 2). 
Equation (28) can also be written as (dropping the 
subscripts) 


and the measurement matrix, H rar<? b as 

^HfAdTaxPixsl (24) 40 

results in the following linear update equation 

y<,=H rate , i £ rate , b +T\ (25) 

The innovations for the rate estimation are then given as 45 


z= Q r §j} = Q^Ds J (29) 

where 

a=baseline direction expressed in body coordinates; 
D=direction cosine matrix, transforms from inertial to 
body coordinates; and 

s 7 =satellite to GPS direction expressed in inertial coordi- 
nates. 

Taking the derivative of (29) results in the following 


Z - yio Hrate,blL mte b 


(26) 


Z = a T bs[ + a T DS[ 


(30) 


where 

K ral e b = the estimated rate 

and b state elements. 

The elements in y^, given in (23), are computed numeri- 
cally. The derivative of the line of sight, 

s, 

can be computed from the given GPS position and velocity 
and the estimated position and velocity. The angle, 0, is 
computed from the phase measurement and is numerically 
differentiated. This defines the derivative approach for GPS. 65 
Alternatively, the ‘estimation approach’ could be utilized for 
the GPS rate estimation. 


50 The time derivative of the direction cosine matrix can be 
written as 

Z)=-[to^]D (31) 

55 Where 

[cox]=cross product matrix of the rotation rate. 

Substituting (31) into (30) and reordering the vectors 
gives 

60 z = a T [s b <r-]u) + a T bs, ^ 

where 

bfc x ]=[Ds 7 x]=cross product matrix of the GPS direction 
unit vector in body coordinates; 

D=current estimate of the direction cosine matrix which 
transforms from inertial to body coordinates; and 
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s 7 =derivative of the GPS direction unit vector in inertial 
coordinates, this is not computed numerically but is 
developed from the relative position and velocity vec- 
tors. 

Development of the GPS Line of Sight Derivative is 5 
discussed here. The position vector from the spacecraft to a 
GPS satellite is given as 


z = F M + u (39) 

Equation (39) is used to propagate z numerically. 

In general, the linearized error dynamics are given as 

x{t) = Fx(t) (40) 


R=R SV -R 


where 


(33) 

10 


R sv =the position vector of the GPS satellite from the 
center of Earth; and 

R=the position vector of the satellite from the center of 15 
Earth. 

The GPS line of sight unit vector is then given as 


R s 

i/£i 


(34) 

20 


where 

F=ff//X 

The full dynamics matrix, F, is given above as well as in 
Reference 1 for all the states except the phase measurement 
states, z. The matrix Fz from (39) into augmented into the 
matrix F as 

F orbit 0 0 0 - 

_ 0 F clock 0 0 

0 0 Fattitude 0 

.0 0 0 F rateAz _ 


Taking the derivative of s results in 


d 1 V/ — 1 


( 35 ) 25 


The filter covariance matrix is propagated according to 


flk + i(-) = «i(i*(+))P* (+)*%(+)) + 0t 


(41) 


The first term on the right hand side of (35) can be written 
as 


dr i R s -v s (36) 

rfrliiyj"” I Rs\ 3 


where 

V^=V^ v -V=R^=the relative velocity of the GPS satellite 
with, respect to the spacecraft 
Substituting (36) into (35) results in 


; R S 7V S 1 

S ~~~iRF Rs + \RA Vs 


(37) 


The GPS position and velocity vectors R sv and V sv , 
respectively, are provided. However, only estimates of the 
spacecraft position and velocity vectors R and V, 
respectively, are available. Resolving all the vectors in 
inertial coordinates, and substituting in the estimated values 
results in the following expression for the derivative of the 
GPS line of sight vector implemented into the EKF algo- 
rithm. 


where 

30 <D=the state transition matrix=e jFAr ; and 

Q=matrix of standard deviations of the white noise pro- 
cesses driving each of the states. 

The filter uses the ‘multiplicative’ approach in the attitude 
estimation 3 , therefore the error state for the attitude is a 
35 vector of attitude errors. 

Filter Update 

The measurement matrix, H, is partitioned as \\=\\\ orbit , 
H attitude 9 Hrate , b The individual measurement matrices are 

defined according to the two measurements used, namely the 
40 magnetometer measurement and the GPS phase and range 
measurements. The development of the magnetometer and 
GPS measurement matrices are above. The development of 
the GPS measurement matrix for the position and velocity 
states as above includes both range and range-rate type 
45 measurements. In the current version of the EKF, only range 
measurements are included. 

GPS Phase Measurement State Update 

According to the development of the ‘estimation 
approach’ given in Reference 4, the measured phase z m9 is 
50 related trivially to the rate and the state z as 


= [0 0 /] 


■*" XL ~ ^ rate,b,z^rate,b,z 4" 2 


(42) 


R s iVsi- 1 * 

-^L^r sJ + -—v sJ 
I'M \R s ,i I 


(38) 


where the "symbol refers to estimated values, and the ‘I’ 
refers to inertial coordinates. 


where 

r|.=added white noise 

The innovations are then computed as 

Z.=Zm-H rate&’ZKrcite&z 


(43) 


The vector z is formed from the individual phase 
measurements, z tj (either two components for the two base- 
lines and one GPS satellite or four components for the two 
baselines and two GPS satellites, or two for each GPS 
satellite observed if more than two are observed). Equation 
(32) is then written in the general form as 


and are used in the update equations for a linear Kalman 
filter. The ^ rate byZ matrix replaces the portion of the GPS 
measurement matrix developed for the ‘derivative approach’ 
65 given in Reference 5. 

The use of the ‘estimation approach’ for GPS requires the 
initialization of the z states whenever a new GPS satellite is 
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being observed. The states are initialized a priori with the 
first phase observations. The covariance is initialized with 
the a priori values and the correlation terms are set to zero. 
Simulations and Results 
Case 1 

The test scenario consisted of a simulated orbit for a low 
earth orbit spacecraft and simulated orbits for the GPS 
constellation. The spacecraft rotated at a 1 revolution per 
orbit rate about the y body axis in an earth pointing attitude, 
at an altitude of 700 km and inclination of 98 degrees. The 
two GPS antennas are aligned along the spacecraft x and y 
body axes, on the anti-nadir side of the spacecraft. The 
magnetometer is aligned along the spacecraft body axes. The 
magnetometer and reaction wheel data contained an added 
white noise, the GPS measurements had no corruption. All 
the data was output at a rate of 1 Hz. The errors between the 
true state and the a priori filter state are given in Table I, 
except for the state, b. The state, b, was initialized with the 
first magnetometer measurement. The true state of the filter 
is output at 1 Hz also. 


TABLE 1 


A Priori Errors 

STATE 

in the Filter State 
ERROR 

Position 

100 km/axis 

Velocity 

0.5 km/sec/axis 

Quaternion 

12 deg rotation error 

Rate 

Orbit rate in the pitch 


axis only 


Results 

The initial test of the GPS/magnetometer navigation sys- 
tem included 30 minutes of simulated data. FIGS. 2 through 
5 show the results of this first test. All the state parameters 
plotted appear to be converging after approximately 10 
minutes. All of the states experience large deviations during 
the first 10 minutes. FIG. 2 shows the position errors about 
each axis. By 30 minutes the errors are approximately 1 km 
in each axis. FIG. 3 shows the velocity errors about each 
axis. FIG. 4 shows the angular attitude error about each of 
the body axes. By 30 minutes the errors are reduced to less 
than 1 degree in each axis. Finally, FIG. 5 shows the rate 
errors about each axis. The rate errors in the x and y body 
axes look noisier than those in the z body axis. A plot of the 
additional state that was added according to the estimation 
approach of is not included. It follows the magnetic field 
vector in body coordinates, as expected. In this preliminary 
test case, the filter included data from the magnetometer and 
only one GPS satellite per update cycle. Additional GPS 
satellites in view could easily be processed per update cycle. 
Advantages of processing only one per cycle are the reduc- 
tion in processing time per update and in the reduction of 
data that must be carried between two update cycles (for use 
in the derivative computations). 

Case 2 

The test scenario consisted of an orbital ephemeris of the 
Upper Atmosphere Research Satellite (UARS) spacecraft 
and simulated orbits for the GPS constellation. The space- 
craft is inertially pointing on an orbit whose average altitude 
is 575 km, its inclination is 57 degrees, and its eccentricity 
is 0.0015. The simulation generated from the UARS ephem- 
eris consists of 12 hours of sensor and reaction wheel data. 
After 6 hours, an attitude maneuver is simulated about the y 
axis at a rate of 0.2 deg/sec for 6.63 minutes (giving a total 
maneuver of 80 degrees). Two GPS antennas are aligned 
along the spacecraft x and y body axes and one at the origin, 
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pointing in the -z body axis direction. A magnetometer is 
aligned along the spacecraft body axes. The simulated 
magnetometer and GPS phase and range data contain added 
white noise, with standard deviations of 1 milliGauss, 1 
5 degree, and 2 m, respectively. All the data are output at a rate 
of 2 Hz. The initial errors between the true state and the a 
priori filter state are given in Table I, except for the states b 
and z. The states b and z are initialized with the first 
magnetometer and GPS phase measurements, respectively. 
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TABLE 2 


INITIAL A PRIORI ERRORS IN THE FILTER STATE 

STATE 

ERROR 

Position 

1100 km (RSS) 

Velocity 

0.6 km (RSS) 

Quaternion 

111 deg rotation error 

Rate 

0.06 deg/sec/axis 


20 Results 

FIGS. 6 and 7 show the RSS of the angular velocity errors 
for measurements from one and two GPS satellites, 
respectively, during the first 150 minutes of the data span. In 
both cases the estimation approach is used (as applied to the 
25 GPS phase measurements, recalling that the estimation 
approach is applied to the magnetometer). Here it can be 
seen that the use of two GPS satellites shortens the conver- 
gence time. After the transients have settled there is little 
difference between the two cases. The RSS rate errors using 
30 the derivative approach are slightly noisier, but show similar 
behavior to the estimation approach for both the one GPS 
and two GPS satellite scenarios. The jumps in the RSS errors 
within the first 50 minutes occur when there is a change to 
a different GPS satellite. 

35 FIGS. 8 and 9 show the estimated y angular velocity for 
the estimation and the derivative approaches, respectively, 
each using measurements from one GPS satellite. The time 
span in the Figures is during the attitude maneuver, which 
occurs after 6 hours. While in FIGS. 6 and 7 we compared 
40 the effect of one and two GPS satellites (because the 
difference between the estimation and derivative approaches 
was negligible), in FIGS. 8 and 9 we compare the two 
approaches because here the difference stems from the 
approaches and the difference between using one or two 
45 GPS satellites is negligible. In both cases, the filter follows 
the maneuver. The other two axes, namely the x and z axes, 
stay approximately at 0 deg/sec, similar to the estimated y 
rate both before and after the maneuver. 

FIGS. 10 and 11 show the RSS attitude error during the 
50 first 150 minutes, using measurements from one and two 
GPS satellites, respectively, using the estimation approach. 
Again, the use of two GPS satellites shortens the conver- 
gence time. After the filter has converged the results are 
similar with final average errors of approximately 0.3 
55 degrees (computed over the last 30 minutes). The derivative 
approach results are similar to the estimation approach 
results. 

Finally, FIGS. 12 and 13 show the RSS position errors for 
measurements from one and two GPS satellites, 
60 respectively, using the estimation approach. Here, the dif- 
ference between using one or two GPS satellites is the most 
significant (the same is true for the velocity errors which 
have plots similar in shape to the position errors). After 12 
hours, the one GPS case is still converging, with peak errors 
65 of about 4.5 km and an average RSS over the last 30 minutes 
of 0.8 km. The spikes in the errors occur due to a change in 
the GPS satellite and there is also an error from the magnetic 
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field that occurs at the orbital period (approximately 90 
minutes). In the case of using measurements from two GPS 
satellites, the effects of the changing GPS satellites are 
reduced. After 12 hours, the peak errors are reduced to less 
than 0.3 km with an average over the last 30 minutes of 30 5 
m. The results for the final average RSS velocity errors are 
6.6 cm/sec using two GPS versus 100 cm/sec when using 
one GPS. 

What is claimed is: 

1. An autonomous navigation system for determining a set 
of satellite navigation parameters including attitude, orbit, 
and rate, the system comprising: 

a GPS receiver; 

a magnetometer; 

processing means for processing measurements from the 
GPS and magnetometer to estimate attitude, orbit, and 
rate. 

2. The autonomous navigation system as in claim 1, 
wherein the measurements from the GPS comprises: 

GPS pseudo -range; and GPS phase measurements. 

3. The autonomous navigation system as in claim 2, 20 
wherein the measurements from the GPS further comprises 
GPS range rate. 

4. The autonomous navigation system as in claim 1, the 
processing means comprising: 

an algorithm comprising a Kalman filter. 

5. The autonomous navigation system as in claim 4, 
wherein output filter state of the Kalman filter includes: 

x T =[R T , v T , c T , a T , 0 /], 

30 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, g is the attitude quaternion, co is the 
rotation rate. 

6 . The autonomous navigation system as in claim 5, 35 
wherein the rate estimate is obtained by applying derivative 
approach to GPS phase measurements or magnetometer 
measurements. 

7. The autonomous navigation system as in claim 5, 
wherein the rate estimate is obtained by applying derivative 40 
approach to GPS phase measurements and magnetometer 
measurements. 

8 . The autonomous navigation system as in claim 4, 
wherein output filter state of the Kalman filter includes: 

45 

X t =[R t , V t , C t , a T , co r , b T l 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, g is the attitude quaternion, co is the 50 
rotation rate, and b is the magnetic field in body 
coordinates. 

9. The autonomous navigation system as in claim 8 , 

wherein the rate estimate is obtained by applying derivative 
approach to GPS phase measurements. 55 

10. The autonomous navigation system as in claim 4, 
wherein output filter state of the Kalman filter includes: 

x t =[r t , v T , c T , a T , co r , f], 

60 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, g is the attitude quaternion, co is the 
rotation rate, and z is GPS phase measurement. 

11. The autonomous navigation system as in claim 10, 65 
wherein the rate estimate is obtained by applying derivative 
approach to magnetometer measurements. 
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12. The autonomous navigation system as in claim 4, 
wherein output filter state of the Kalman filter includes: 

x t =[r t , v t , c t , g T , to 7 , b T , f], 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, b is the attitude quaternion, co is the 
rotation rate, b is the magnetic field in body 
coordinates, and z is GPS phase measurement. 

13. A method for estimating a set of navigation parameters 
of a satellite including attitude, orbit, and rate, the method 
comprising: 

providing measurements from GPS; 

providing measurements from magnetometer; 

executing an algorithm to estimate attitude, orbit, and rate 
of the satellite based on the measurements from the 
GPS and magnetometer. 

14. The method as in claim 13, wherein the measurements 
from the GPS comprises: 

GPS pseudo-range; and 

GPS phase measurements. 

15. The method as in claim 14, wherein the measurements 
from the GPS further comprises GPS range rate. 

16. The method as in claim 13, wherein the step of 
executing the algorithm includes performing Kalman filter- 
ing. 

17. The method as in claim 16, wherein output filter state 
of the Kalman filter includes: 

x t =[r t , v T , c T , a T , to 7 ], 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, g is the attitude quaternion, co is the 
rotation rate. 

18. The method as in claim 17, wherein the rate estimate 
is obtained by applying derivative approach to GPS phase 
measurements or magnetometer measurements. 

19. The method as in claim 17, wherein the rate estimate 
is obtained by applying derivative approach to GPS phase 
measurements and magnetometer measurements. 

20. The method as in claim 16, wherein output filter state 
of the Kalman filter includes: 

x t =[r t , v t , c t , a T , < o r , b 7 ], 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors, g is the attitude quaternion, co is the 
rotation rate, and b is the magnetic field in body 
coordinates. 

21. The method as in claim 20, wherein the rate estimate 
is obtained by applying derivative approach to GPS phase 
measurements. 

22. The method as in claim 16, wherein output filter state 
of the Kalman filter includes: 

x t =[r t , v t , c t , a T , m r , f\ 

wherein R and V are the spacecraft position and velocity 
vectors, respectively, C is a vector of GPS receiver 
clock errors g is the attitude quaternion, co is the 
rotation rate, and z is GPS phase measurement. 
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